% Dynamic bank run model 
% Bank leverage choice only
% GHH preferences + TFP shock

%-------------------------------------------------------------------------------
% 0. Housekeeping, paths, and model
%-------------------------------------------------------------------------------
%set(0,'DefaultFigureWindowStyle','docked')
dbstop if error
dbclear if error

%-------------------------------------------------------------------------------
% 1. Declaration of variables and parameters
%-------------------------------------------------------------------------------
%var cU,yU,iU,kU,hU,nU,LU,RbU,rkhatsU,aU,xiU,PU; % 12 variables
var     yU,iU,xiU,PU,psi_p,n0_p,gam_p,kU,nU,LU,RbU,rkhatsU,aU,cU,hU; % 12 variables
%var psi_p,n0_p,gam_p;                           % 3 parameters to be calibrated

varexo ea en;                                

parameters hU_SS,LU_SS,PU_SS;                 % 3 target values
parameters beta_p,nu_p,alpha_p,delta_p,lam_p,chi0_p,chi1_p,rhoa_p,siga_p; % standard parameters 

%--------------------------------------------------------------------------------
% 2. Read parameters
%--------------------------------------------------------------------------------
@# include "params.mod"

%--------------------------------------------------------------------------------
% 3. Model
%--------------------------------------------------------------------------------
model;
% Economic relationships
%---------------------------------------------------------------------------------
  % Eq.1: Crisis probability
    # murkU = log(alpha_p*((1-alpha_p)/psi_p)^((1-alpha_p)/(alpha_p+1/nu_p)))+(1+nu_p)/(alpha_p*nu_p+1)*rhoa_p*aU-(1-alpha_p)/(alpha_p*nu_p+1)*log(kU);
    # sigrk_p = (1+nu_p)/(alpha_p*nu_p+1)*siga_p;
    PU = (1/2)*(1+erf((rkhatsU-murkU)/(sigrk_p*sqrt(2))));
  
  % Eq.2: Household, Euler equation
    # pi = 3.141592653589793;
    # ERkp1d = exp(murkU+sigrk_p^2/2)*1/2*(1+erf((rkhatsU-murkU-sigrk_p^2)/(sigrk_p*sqrt(2))))+(1-delta_p)*PU;
    1/(beta_p*RbU) = ((cU-psi_p*hU^(1+1/nu_p)/(1+1/nu_p))/(cU(+1)-psi_p*hU(+1)^(1+1/nu_p)/(1+1/nu_p)))*(ERkp1d*LU/(RbU*(LU-1))-lam_p*PU+1-PU);

  % Eq.3: Labor market equilibrium condition
    yU = psi_p/(1-alpha_p)*hU^(1/nu_p+1);    

  % Eq.4. Production   
    yU = exp(aU)*(kU(-1))^(alpha_p)*hU^(1-alpha_p);

  % Eq.5: Law of motion for capital
    kU = (1-delta_p)*kU(-1) + iU;    

  % Eq.6: Good market clearing
    yU = cU + iU;    

  % Eq.7: Bank balance sheet
    kU = LU*nU;    

  % Eq.8: Thresdhold
    exp(rkhatsU) = RbU*(1-1/LU)*(1+lam_p*(1-gam_p))-1+delta_p;    

  % Eq.9: Bank leverage
    # ERkp1u = exp(murkU+sigrk_p^2/2)*(1/2)*(1+erf((murkU+sigrk_p^2-rkhatsU)/(sigrk_p*sqrt(2))))+(1-delta_p)*(1-PU);
    RbU*(1-PU) = ERkp1u - lam_p*RbU^(2)*(1-gam_p)*(1+lam_p*(1-gam_p))*(LU-1)/(LU^2) / (RbU*(1-1/LU)*(1+lam_p*(1-gam_p))-1+delta_p)*1/(sqrt(2*pi)*sigrk_p)*exp(-1/2*((rkhatsU-murkU)/sigrk_p)^2);

  % Eq.10: Law of motion for bank net worth
    nU = chi1_p*(chi0_p*((alpha_p*yU/kU(-1)+1-delta_p)*kU(-1)-RbU(-1)*(1-1/LU(-1))*kU(-1)-nU(-1)) + nU(-1)) + n0_p - en;    

  % Eq.11: TFP shock
   aU = rhoa_p*aU(-1) + siga_p*ea;  

  % Eq.12: Transformed shock
  xiU = alpha_p*yU/kU(-1)/exp(rkhatsU(-1));

  % Calibrated parameters
  psi_p = steady_state(psi_p);
  n0_p = steady_state(n0_p);
  gam_p = steady_state(gam_p);
  
end;


%-----------------------------------------------------------------------------
% 4. Simulations
%-----------------------------------------------------------------------------
steady;

check;

shocks;
var ea;      stderr 1;
var en;      stderr 0.6;
end;

stoch_simul(order = 1,nograph,irf=80); %, periods=200, noprint, nofunctions